% Scrpits to compare LS and ROBUST estimate of amplitudes
% The plan is also to compare pmtm with the mean square spectrum 

% generate a sinusoid with the following parametrs

numsamples = 1000;
Frq = [1/25 1/14]; %frequencies
A = [2 2 5]; %amplitudes
fs = 1/5;%~ 5 minutes of Julia Vz

x = (1:numsamples)/fs;
T = (1/fs) * length(x); % size of the fft bin
y = A(1)*sin(2*pi*x*Frq(1)) + A(2)*sin(2*pi*x*Frq(2) + pi/4) + normrnd(0,A(3),1,length(x));

% introduce outliers on 30 % of the data by randomly adding spikes

k = floor(rand([1,floor(numsamples*0.2)])*numsamples) + 1;

y(k) = normrnd(0,A(3)*20,1,length(k));


% calculate mean-power at the above harmonics by LS fitting

% model
m = [cos(2*pi*x*Frq(1))'       sin(2*pi*x*Frq(1))' ... 
     cos(2*pi*x*Frq(2))'       sin(2*pi*x*Frq(2))'];
 
%ls fitting

spec = m\y'; 

% estimating amplitudes
amp1 = sqrt(spec(1)^2+spec(2)^2);
amp2 = sqrt(spec(3)^2+spec(4)^2);

%The estimated amplitudes amp1 and amp2 should be very close to A(1) and
%A(2) resp.

% robustfitting

spec = robustfit(m,y');

amp3 = sqrt(spec(2)^2+spec(3)^2);
amp4 = sqrt(spec(4)^2+spec(5)^2);

fprintf('LS  expected %d %d estim %f %f\n', A(1),A(2),amp1,amp2);
fprintf('ROB expected %d %d estim %f %f\n', A(1),A(2),amp3,amp4);

% In all the cases, robust fit consistently gave better estimates than LS
% Evene when the % of outliers (defined as random data > 20 times the
% amplitude of the background noise". When there are no "ouliers present in
% the data", LS and ROB are comparable.
% Take home message. Always use robust processing


% Is there a difference when inverting one component and when inverting
% with all the components ?

% two components

% model
m = [cos(2*pi*x*Frq(1))'       sin(2*pi*x*Frq(1))' ... 
     cos(2*pi*x*Frq(2))'       sin(2*pi*x*Frq(2))'];
spec = robustfit(m,y');

amp3 = sqrt(spec(2)^2+spec(3)^2);
amp4 = sqrt(spec(4)^2+spec(5)^2);


% only one component

% model
m = [cos(2*pi*x*Frq(1))'       sin(2*pi*x*Frq(1))'];

spec = robustfit(m,y');

amp5 = sqrt(spec(2)^2+spec(3)^2);

m = [cos(2*pi*x*Frq(2))'       sin(2*pi*x*Frq(2))'];

spec = robustfit(m,y');

amp6 = sqrt(spec(2)^2+spec(3)^2);

fprintf('ROB together %d %d estim %f %f\n', A(1),A(2),amp3,amp4);
fprintf('ROB separate %d %d estim %f %f\n', A(1),A(2),amp5,amp6);

% Inference. The amplitude estimates were closer to the original when
% ALL the two components were included. The difference is more noticeable
% when the noise is high


% I also tested pmtm. But could not get a way to estimate the amplitude
% from the power spectral density